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Abstract 

We present a novel approach for the derivation of PDEs modeling curvature- 
driven flows for matrix-valued data. This approach is based on the Riemannian 
geometry of the manifold of Symmetric Positive Definite Matrices Vin). The 
differential geometric attributes of 'P(n) —such as the bi-invariant metric, the co- 
variant derivative and the Christoffel symbols— allow us to extend scalar- valued 
mean curvature and snakes methods to the tensor data setting. Since the data live 
on 'P{n), these methods have the natural property of preserving positive definite- 
ness of the initial data. Experiments on three-dimensional real DT-MRI data show 
that the proposed methods are highly robust. 

1 Introduction 

With the introduction of diffusion tensor magnetic resonance imaging (DT-MRI) ID, 
there has been an ever increasing demand on rigorous, reliable and robust methods 
for the processing of tensor-valued data such as the estimation, filtering, regularization 
and segmentation. Many well established PDE-based methods used for the processing 
of scalar-valued data have been extended in various ways to the processing of multi- 
valued data such as vector-valued data and smoothly constrained data Il5l fmi20l 12111221 
l23l . Recently, some efforts have been directed toward the extension of these methods 
to tensor fields 13] [8] |7] |T3] [161 HH HS). The generalization of the methods used for 
scalar- and vector-valued data to tensor-valued data is being pursued with mainly three 
formalisms: the use of geometric invariants of tensors like eigenvalues, determinant, 
trace; the generalization of Di Zenzo's concept of a structure tensor for vector-valued 
images to tensor-valued data; and recently, differential-geometric methods. 

The aim of the present paper is to generalize the total variation (TV) flow, mean cur- 
vature motion (MCM), modified mean curvature flow and self snakes to tensor-valued 
data such as DT-MRI. The key ingredient for these generalizations is the use of the 
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Riemannian geometry of the space of symmetric positive-definite (SPD) matrices. The 
remainder of this paper is organized as follows. In Section [2] we give a compilation 
of results that gives the differential geometry of the Riemannian manifold of symmet- 
ric positive-definite matrices. In Section[3]we fix notation and recall some facts about 
immersions between Riemannian manifolds and their mean curvature. We explain in 
Section |4] how to describe a DT-MR image by differential-geometric concepts. Sec- 
tion |5] is the key of our paper in which we extend several mean curvature-based flows 
for the denoising and segmentation from the scalar and vector setting to the tensor one. 
In Section[6]we present some numerical results. 

2 Differential Geometry of V{n) 

Positive-definite matrices are omnipresent in many engineering and physical contexts. 
They play important roles in various disciplines such as control theory, continuum me- 
chanics, numerical analysis, covariance analysis, signal processing, etc. Recently, they 
gained an increasing attention within the diffusion tensor magnetic resonance imag- 
ing (DT-MRJ) community as they are used as an encoding for the principal diffusion 
directions ans strengths in biological tissues. 

We here recall some differential-geometric facts about the space of symmetric 
positive-definite matrices that have been recently published by the authors. We de- 
note by S{ti) the vector space of n x n symmetric matrices. A matrix A E S{n) is 
said to be positive semidefinite if x^Ax > for all x g K", and positive definite if in 
addition A is invertible. The space of all n x n symmetric, positive-definite matrices 
will be denoted by 'P{n). We note that the set of positive-semidefinite matrices is a 
pointed convex cone in the linear space of n x rt matrices, and that is the interior 
of this cone. It is a differentiable manifold endowed with a Riemannian structure. The 
tangent space to 1^(71) at any of its points P is the space TpV{n) ~ {P} x S{n), which 
for simplicity is identified with S{n). On each tangent space TpV{n) we introduce the 
base point-dependent inner product defined by {A, B) p := tr{P~^AP^^B). 

This inner product leads to a natural Riemannian metric on the manifold V{n) that 
is given at each P by the differential 

ds^ = tr {P- 1 dPP- ^dP) , (1) 

where dP is the symmetric matrix with elements (dPij). We note that the metric 
([1} is invariant under congruent transformations: P LPL^ and under inversion 

P~^ p-\ 

For an 71 X n matrix A we denote by vec A the n^-column vector that is obtained by 
stacking the columns of A. If A is symmetric, then then ^n{n — 1) elements of vec{A) 
are redundant. We will denote by v{A) the d — |n(n+ l)-vector that is obtained from 
vec{A) by eliminating the redundant elements, e.g., all supradiagonal elements of A. 
We note that there are several ways to arrange the independent elements of vec (A) into 
v{A). In any case, there exists a unique x ^n{n + 1) matrix, called the duplication 
matrix and denoted by £)„, that by duplicating certain elements, reconstructs vec A 
from v{A), i.e., is the matrix such that 

vecA^ D„v{A). (2) 
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The duplication matrix which has been studied extensively by Henderson and 
Searle |9|, and by Magnus and Neudecker [141, has full column rank ^n{n+l). Hence, 
Dj^Dn is non-singular and it follows that the duplication matrix £)„ has a Moore- 
Penrose inverse denoted by _D+ and is given by 

Dl^{DlDr)^' Dl. 

It follows from (|2| that 

v{A)^D+vecA. (3) 

By using the vector v{P) as a parametrization of P e V{n) we obtain the matrix 
of components of the metric tensor associated with the Riemannian metric ([T]i is given 
explicitly by 1.25 J 

G{P) = DI{P-^ ® P-^)D^. (4) 

For differential-geometric operators on V{n) it is important to obtain the expression of 
the inverse of the metric and that of its determinant. The matrix of components of the 
inverse metric tensor is given by 

G-\P) = D+ {P ® P) of , (5) 

and the determinant of G is 

det(G(P)) = 2"(""i)/2 ((det(P))("+i) . (6) 
In the coordinate system (p"), the Christoffel symbols are given by 1251 

= {P~' ® E^) Dr^U, 1 < a, /3, 7 < rf, 

where E"^ is the dual basis associated with the local coordinates {p"). As the elements 
of E"' and Z)„ are either 0, 1, or |, it follows from the above theorem that each non- 
vanishing Christoffel symbol is given by an element of P^^ or half of it. 

Let P be an element of 7^(3) and let dP be a (symmetric) infinitesimal variation of 
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Hence, the complete and reduced vector forms of P are respectively, 

vec(P) = y / / / p2 pb p6 p5 p3]T^ ^(p) ^ ^pl p2 p3 p4 ^5 ^ejT^ 

The components of the inverse metric tensor and the Christoffel symbols are given 
explicitly in the appendix. 
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3 Immersions and Mean Curvature 



Let (M, 7) and {N, g) be two connected Riemanian manifolds of dimensions m and n, 
respectively. We consider a map (f) : M ^ N that is of class C^, i.e., cj) e C^(Af, A^). 
Let {x"} — 1 < a < mhea local coordinate system of a; in a neighborhood of a point 
p E M and let {y'}i<,;<ri be a local coordinate system of y in a neighborhood of 
(j){P) EN. 

The mapping induces a metric 0*5 on Ad defined by 

c^*g{Xj,,Yp)^g{cl>,{Xj,),cl>,{Yp)). (7) 

This metric is called the pull-back metric induced by as it maps the metric in the 
opposite direction of the mapping 0. 

An isometry is a diffeomorphism (f> : M ^ N that preserves the Riemannian 
metric, i.e., if g and 7 are the metrics for M and N, respectively, then 7 = (j>*g. It 
follows that an isometry preserves the length of curves, i.e., if c is a smooth curve on 
M, then the curve o c is a curve of the same length on N. Also, the image of a 
geodesic under an isometry is again a geodesic. 

A mapping : M — > iV is called an immersion if (0*)^ is injective for every point 
p in M. We say that M is immersed in by or that M is an immersed submanifold 
of N. When an immersion (p is injective, it is called an embedding of M into N. We 
then say that M is an embedded submanifold, or simply, a submanifold of N. 

Now letcf) : M ^ N be an immersion of a manifold M into a Riemannian manifold 
N with metric g. The first fundamental form associated with the immersion (p is h = 
(j)*g. Its components are haf3 = daCp^dpcj)^ gij where daij)^ — The total covariant 
derivative Vd(j) is called the second fundamental form of (p and is denoted by 11^^ {</)). 
The second fundamental form 11^^ takes values in the normal bundle of M. The mean 
curvature vector H of an isometric immersion (p : M ^ N is defined as the trace of 
the second fundamental form 11^^ {(f>) divided by m = dim AI ifTOl 

H —tr^II^Ucj)). (8) 
m 

In local coordinates, we have ifTOl 

= Am0^ + r^ix)-T), W-)) 1^1^. (9) 

where ^r* j, are the Christoffel symbols of {N, g) and Am is the Laplace-Beltrami 
operator on (M, 7) given by 



4 Diffusion-Tensor MRI Data as Isometric Immersions 

A volumetric tensor-valued image can be described mathematically as an isometric 
immersion {x^^x^^x^) = (a;^, x^, x"^; P(a;^, x^, a;'^)) of a three-dimensional do- 
main Q. in the fiber bundle M"^ (g) 7^(3), which is a nine-dimensional manifold. We 
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denote by {Ad, 7) the image manifold and its metric and by (iV, g) the target manifold 
and its metric. Here Af = 17 and N — ® ^(3)- Consequently, a tensor-valued 
image is a section of this fiber bundle. The metric g of is given by 

ds^ = dSgpatial + '^^tcnsoi- (H) 

The target manifold N, in this context is also called the space-feature manifold 
II23I . We can rewrite the metric defined by ( [TT| as the quadratic form 

df = {dx^f + {dx'Y + {dx^f + {dpfDliP-^ ® p-^)D„{dp), 

where p — (p^) — v{P). The corresponding metric tensor is 

^3 03,6 

06,3 9 

where g is the metric tensor of 7^(3) as defined in Section |2] 

Since the image is an isometric immersion, we have — (j>*g. Therefore 

lap^ 5aP+ gijdaP'dpp' , tt, /3 = 1, . . . , TO, = (12) 

We note that d = n — mis the codimension of M. In compact form, we have 

7 - /„, + iVpfGi^)Vp. (13) 

where G is given by Q. (We take to = 2 for a slice and to = 3 for a volumetric 
DT-MRI image.) 



5 Geometric Curvature-Driven Flows for Tensor- Valued 
Data 

The basic concept in which geometric curvature-driven flows are based is the mean 
curvature of a submanifold embedded in a higher dimensional manifold. Here we 
generalize the scalar mean curvature flow to mean curvature flow in the space-feature 
manifold (g) P(3). For this, we embed the Euclidean image space into the Rie- 
mannian manifold il (g) 7^(3), and use some classical results from differential geometry 
to derive the Riemannan Mean Curvature (RMC). We then use the RMC to generalize 
mean curvature flow to the tensor-valued data. Given the expression of the mean cur- 
vature vector H, we can establish some PDEs based tensor-image filtering. Especially, 
we are interested of the so called level-set methods, which relay on PDEs that modify 
the shape of level sets in an image. 

5.1 Riemannian Total Variation Flow 

The total variation norm (TV) method introduced in |I91 and its reconstructions have 
been successfully used in reducing noise and blurs without smearing sharp edges in 
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grey-level, color and other vector- valued images |I6] [12] [TTl |20l IIU |22]| . It is then 
natural to look for the extension of the TV norm to tensor- valued images. 

The TV norm method is obtained as a gradient-decent flow associated with the L^- 
norm of the tensor field. This yields the following PDE that express the motion by the 
mean curvature vector H 

dtcj)' = H\ (14) 

This flow can be considered as a deformation of the tensor field toward minimal im- 
mersion. Indeed, it derives from variational setting that minimize the volume of the 
embedded image manifold in the space-feature manifold. 



5.2 Riemannian Mean Curvature Flow 

The following flow was proposed for the processing of scalar-valued images 

Vu 

dtu=\'^u\d\Y —— u{Q,x,y) = uo{x,y), (15) 
|Vu| 

where uq{x, y) is the grey level of the image to be processed, u{t, x, y) is its smoothed 
version that depends on the scale parameter t. 

The "philosophy" of this flow is that the term | Vu| div represents a degenerate 
diffusion term which diffuses u in the direction orthogonal to its gradient Vw and does 
not diffuse at all in the direction of Vu. 

This formulation has been proposed as a "morphological scale space" fS) and as 
more numerically tractable method of solving total variation ifTSl . 

The natural generalization of this flow to tensor- valued data is 

dtc/)' ^\W^cl>\gH\ i^l,...,d. (16) 

where 

We note that several authors have tried to generahze curvature-driven flows for tensor- 
valued data in different ways. We think that the the use of differential-geometric tools 
and concepts yield the correct generahzation. 



5.3 Modified Riemannian Mean Curvature Flow 

To denoise highly degraded images, Alvarez et al. llil have proposed a modification of 



the mean curvature flow equation ( 15 i that reads 



at0 = c(|if*V0|)|V0|div^, <f>iO,x,y) = Mx,y), (17) 

where if is a smoothing kernel (a Gaussian for example). A' ★ V0 is therefore a local 
estimate of V0 for noise elimination, and c{s) is a nonincreasing real function which 
tends to zero as s oo. We note that for the numerical experiments we have used 

c(|V0|) = fcV(fc' + |V0p). 
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The generalization of the modified mean curvature flow to tensor-field processing 

is 

dt^' ^c{\K^V'^ct>\a)\^^^\aH\ 4>\Q,^)^4>m. (18) 

The role of c is to reduce the magnitude of smoothing near edges. In the scalar case, 
this equation does not have the same action as the Perona-Malik equation of enhancing 
edges. Indeed, Perona-Malik equation has variable diffusivity function and has been 
shown to selectively produce a "negative diffusion" which can increase the contrast 



of edges. Equation of he form (17 1 have always positive or forward diffusion, and 
the term c merely reduces the magnitude of that smoothing. To correct this situation, 
Sapiro have proposed the self-snakes formalism ll20l . which we present in the next 
subsection and generalize to the matrix-valued data setting. 

5.4 Riemannan Self-Snakes 

The method of Sapiro, which he names self-snakes introduces an edge-stopping func- 
tion into mean curvature flow 

= |V0|div(c(if*|V(/.|)^| ^^^^ 
= c{K^ |V0|) I V0I div (^) +Vc{K^ I V0I) • V</) 

Comparing equation (19i to ([iTji, we observe that the term Vc{K -k |V(/)|) • V(f> is 
missing in the old model. This is due to the fact that the Sapiro model takes into 
account the image structure. Indeed, equation ( [T9| can be re-written as 

dt(f> = .?Miffusion + -^^shock, (20) 

where 

^diffusion = c{K*\V(j)\) |V0|div 

^shock = Vc(if*|V0|)-V(^. 

The term .^diffusion is as in the anisotropic flow proposed in |[T]. The second term in 



(20 1, i.e., Vc • V0, increases the attraction of the deforming contour toward the bound- 
ary of "objects" acting as the shock-filter introduced in [18 J for deblurring. Therefore, 
the flow Vc • V0 is a shock filter acting like the backward diffusion in the Perona-Malik 
equation, which is responsible for the edge-enhancing properties of self snakes. See 
i20J for detailed discussion on this topic. 

We are now interested in generalizing Self-Snakes method for the case of tensor- 
valued data. We will start the generaUzation from equation (j20| in the following man- 
ner 

dt4> = T diffusion + ^ shock, (21) 

where 

J'd^ffus^on - C (X * | V^^lg) | V^^lgi?* 

Tshock = Vc(if*|VT</)|g)-VT(/.\ ^^^^ 

This decomposition is not artificial, since the covariant derivative on follow the same 
chain rule as the Euclidean directional derivative: let V a vector field on M which 
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components are w% and let p a scalar function. From the classic differential geometry 
we have 

Vj{pv')^ pVjv' + v'V1p (23) 

and in compact form 

div^ {pV) = p div^ V + V ■ grad^p. (24) 

6 Numerical Experiments 

In Fig. 1 (left), we give a slice of a 3D tensor field defined over a square in M^. We 
note that a symmetric positive-definite 3x3 matrix P is represented graphically by 
an ellipsoid whose principal directions are parallel to the eigenvectors of P and whose 
axes are proportional to the eigenvalues of P^^. Figure 1 (right) shows this tensor field 
after the addition of noise. The resultant tensor field Po{x^ , a;^, a;'^) is used as an initial 
condition for the partial differential equations ( [2T| which we solve by a finite difference 
scheme with Neumann boundary conditions. We used 50 time steps of O.Ol.s. Figure 2 



represents the tensor smoothed by (21 
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Figure 1: Original tensor field (left)and noisy tensor field (right). 



In this paper we generalized several curvature-driven flows of scalar- and vector- 
valued data to tensor-valued data. The use of the differential-geometric tools and 
concepts yields a natural extension of these well-known scalar-valued data processing 
methods to tensor-valued data processing. 
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Appendix 



We give here the explicit form of the inverse metric tensor and Christoffel symbols for 
the Riemannian metric on P{3) . The components of the inverse metric tensor are given 
by 
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The determinant of Pis p = dot P = pVV+2pVp^-P^(p^)^-P^(p^)^-p3(p'^)^r 
Let s := [s^ , s'^ , , s'^ , , s^]"^ = i;(adj(P)), where adj(P) = pP~^ is the adjoint 
matrix of P. 

The Christoffel symbols are arranged in the following six symmetric matrices: 
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